## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.2.1 ✔ dplyr 1.1.2
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
##
## Attaching package: 'lubridate'
##
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
##
##
## corrplot 0.92 loaded
##
##
## Attaching package: 'ggpp'
##
##
## The following object is masked from 'package:ggplot2':
##
## annotate
##
##
##
## Attaching package: 'flextable'
##
##
## The following object is masked from 'package:purrr':
##
## compose
##
##
##
## Attaching package: 'maps'
##
##
## The following object is masked from 'package:purrr':
##
## map
##
##
## ℹ Google's Terms of Service: ]8;;https://mapsplatform.google.com<https://mapsplatform.google.com>]8;;
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
##
## Attaching package: 'scales'
##
##
## The following object is masked from 'package:purrr':
##
## discard
##
##
## The following object is masked from 'package:readr':
##
## col_factor
##
##
## Loading required package: sp
##
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
##
##
## Attaching package: 'raster'
##
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
## rgeos version: 0.6-3, (SVN revision 696)
## GEOS runtime version: 3.10.2-CAPI-1.16.0
## Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## GEOS using OverlayNG
## Linking to sp version: 1.6-0
## Polygon checking: TRUE
##
##
##
## Attaching package: 'rgeos'
##
##
## The following object is masked from 'package:dplyr':
##
## symdiff
This script reads in cleaned datafiles with stream temperature and invertebrate data produced by the invert_data_cleaning.Rmd and teton_data_cleaning.Rmd scripts and performs more involved analyses of stream temperature and invertebrate density/diversity over time. A few findings of interest from previous analyses are that:
* Average summer (July-Aug) stream temperatures have not changed significantly over time
* Average summer (July-Aug) air temperatures have increased over time
* Overall, invetebrate richness has not changed over time
* Invertebrate communities in snow-fed streams appear to be more variable than those in glacier-fed streams
To build on these findings, the main sections of this script are:
I. Temperature:
1. Overall trends in stream temperature over time - Time series analysis: Site-by-site and grouped by source
2. Trends in summer mean, max, min, and degree-days over time: Site-by-site and grouped by source
4. Modeling stream temperature using SNOTEL data: Site-by-site and grouped by source
II. Invertebrates:
1. Trends in richness, abundance, and density over time: Site-by-site and grouped by source
2. Modeling richness, abundance, and density using SNOTEL data: Site-by-site and grouped by source
III. Temperature and Invertebrates
1. Test for correlations between temperature metrics (mean/max/min summer; degree days): Site-by-site and grouped by source
Read in necessary files:
stream_temp <- read.csv("Temperature/cleaned_full_datasets/temps_hourly.csv") #hourly stream temp data
invert_density <- read.csv("invert_data/cleaned_csv/full_invert_densities.csv")%>%
rename(site = Stream)#invert density data
invert_richness <- read.csv("invert_data/cleaned_csv/invert_richness.csv") #invert richness data
snotel <- read.csv("teton_snotel.csv") #snotel data
sources <- read.csv("source_info.csv")%>%
rename(site = stream) #source info, rename for merge
### adding source info to stream temps and invert datasets:
stream_source <- merge(stream_temp, sources, all = T)%>%
mutate(date1 = ymd(date1)) #Add source info to hourly temperature data; re-code date as r date format
stream_source_daily <- stream_source%>%
group_by(site, full_name, stream_code, date1, year, lat, long)%>%
summarise(t_xbar = mean(temp_c, na.rm = T),
t_max = max(temp_c, na.rm = T),
t_min = min(temp_c, na.rm = T),
source = unique(source)) #calculate daily mean, min, and max temps for each site
## Warning: There were 10690 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `t_max = max(temp_c, na.rm = T)`.
## ℹ In group 41: `site = "cloudveil"`, `full_name = "Cloudveil"`, `stream_code =
## "CV"`, `date1 = 2019-09-11`, `year = NA`, `lat = 43.72551`, `long =
## -110.7961`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 10689 remaining warnings.
## `summarise()` has grouped output by 'site', 'full_name', 'stream_code',
## 'date1', 'year', 'lat'. You can override using the `.groups` argument.
density_source <- merge(invert_density, sources, all = T)
richness_source <- merge(invert_richness, sources, all = T)
First, calculate summer (july-aug) mean daily temp, avg daily max and min; plus degree days over 15c? for each site:
stream_summer <- stream_source_daily %>%
mutate(mo = month(date1))%>% #add month column
filter(mo == 7 | mo == 8)%>% #extract July and August
group_by(site, full_name, stream_code, year, lat, long)%>% #group by site and year
summarise(summer_xbar = mean(t_xbar, na.rm = T),
max_xbar = mean(t_max, na.rm = T),
min_xbar = mean(t_min, na.rm= T), #calculate mean temperature, mean daily max and min
dd_2.5 = length(which(t_max >= 2.5)),
dd_5 = length(which(t_max >= 5)),
dd_10 = length(which(t_max >= 10)), #calculate no. days with max temp >2.5, 5, and 10 deg
source = unique(source)) #retain source col
## `summarise()` has grouped output by 'site', 'full_name', 'stream_code', 'year',
## 'lat'. You can override using the `.groups` argument.
Next, extract sites with at least 3 years of data:
under3 <- stream_summer %>%
group_by(site)%>% #group by site
filter(!is.na(year))%>% #remove rows with no year
summarise(n_yrs = length(unique(year)))%>% #count num unique years for site
filter(n_yrs < 3) #extract sites with <3yrs data
summer_3plus <- stream_summer %>%
filter(!site%in%under3$site)
Trends in mean temp:
Nest data by site:
summer_nested <- summer_3plus %>%
filter(!is.na(year))%>%
nest(data = -site)
lm for each site:
summer.mean.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.mean.lms
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 12]> year 0.458 0.309 1.48 0.378
## 2 delta <gropd_df [4 × 12]> year 0.0246 0.0193 1.28 0.330
## 3 grizzly <gropd_df [4 × 12]> year 1.21 0.716 1.68 0.234
## 4 mid_teton <gropd_df [6 × 12]> year 0.0296 0.0280 1.05 0.351
## 5 n_teton <gropd_df [6 × 12]> year 0.255 0.460 0.554 0.609
## 6 paintbrush <gropd_df [5 × 12]> year 1.06 0.510 2.08 0.129
## 7 s_ak_basin <gropd_df [4 × 12]> year 0.0539 0.0159 3.40 0.0768
## 8 s_cascade <gropd_df [3 × 12]> year -0.0705 0.0549 -1.28 0.421
## 9 s_teton <gropd_df [6 × 12]> year 0.987 0.428 2.30 0.0826
## 10 skillet <gropd_df [3 × 12]> year 0.370 0.0617 5.99 0.105
## 11 windcave <gropd_df [3 × 12]> year 0.119 0.0672 1.77 0.328
Nothing significant at 0.05 level; weakly significant increase in mean summer temp at S AK basin (p = 0.08)
Visualization:
summer_mean_withps <- merge(summer_3plus, summer.mean.lms)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
summer.mean <- ggplot(summer_mean_withps, aes(x = year, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.mean
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 row(s) containing missing values (geom_path).
Mean daily max temp
summer.max.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(max_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.max.lms
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 12]> year 0.532 0.328 1.62 0.352
## 2 delta <gropd_df [4 × 12]> year 0.0654 0.0816 0.801 0.507
## 3 grizzly <gropd_df [4 × 12]> year 1.66 1.10 1.51 0.271
## 4 mid_teton <gropd_df [6 × 12]> year -0.0257 0.0485 -0.529 0.625
## 5 n_teton <gropd_df [6 × 12]> year 0.228 0.603 0.378 0.724
## 6 paintbrush <gropd_df [5 × 12]> year 1.74 0.853 2.04 0.135
## 7 s_ak_basin <gropd_df [4 × 12]> year 0.0774 0.108 0.713 0.550
## 8 s_cascade <gropd_df [3 × 12]> year -0.130 0.133 -0.981 0.506
## 9 s_teton <gropd_df [6 × 12]> year 1.37 0.640 2.13 0.0999
## 10 skillet <gropd_df [3 × 12]> year 0.365 0.0254 14.4 0.0442
## 11 windcave <gropd_df [3 × 12]> year 0.158 0.101 1.56 0.364
Significant increase at Skillet, non-significant increases at all but S cascade and mid teton
Visualization:
summer_max_withps <- merge(summer_3plus, summer.max.lms)
summer.max <- ggplot(summer_max_withps, aes(x = year, y = max_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily max. July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.max
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 row(s) containing missing values (geom_path).
Mean daily min temp
summer.min.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(min_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.min.lms
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 12]> year 0.400 0.273 1.46 0.382
## 2 delta <gropd_df [4 × 12]> year 0.0156 0.0152 1.02 0.413
## 3 grizzly <gropd_df [4 × 12]> year 0.849 0.493 1.72 0.227
## 4 mid_teton <gropd_df [6 × 12]> year 0.0639 0.0261 2.45 0.0708
## 5 n_teton <gropd_df [6 × 12]> year 0.254 0.334 0.761 0.489
## 6 paintbrush <gropd_df [5 × 12]> year 0.662 0.325 2.04 0.135
## 7 s_ak_basin <gropd_df [4 × 12]> year 0.0398 0.0440 0.906 0.461
## 8 s_cascade <gropd_df [3 × 12]> year -0.0249 0.0117 -2.14 0.279
## 9 s_teton <gropd_df [6 × 12]> year 0.667 0.292 2.28 0.0846
## 10 skillet <gropd_df [3 × 12]> year 0.347 0.0963 3.61 0.172
## 11 windcave <gropd_df [3 × 12]> year 0.0993 0.0354 2.80 0.218
Non-significant increases at all sites except south cascade
Visualization:
summer_min_withps <- merge(summer_3plus, summer.min.lms)
summer.min <- ggplot(summer_min_withps, aes(x = year, y = min_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily min. July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.min
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 row(s) containing missing values (geom_path).
Degree Days
Days over 2.5 C:
summer.2.5.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(dd_2.5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `model = purrr::map(data, ~lm(dd_2.5 ~ year, data = .) %>%
## tidy)`.
## ℹ In group 8: `site = "s_cascade"`.
## Caused by warning in `summary.lm()`:
## ! essentially perfect fit: summary may be unreliable
summer.2.5.lms
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 12]> year 6.00e+ 0 6.35e+ 0 0.945 0.518
## 2 delta <gropd_df [4 × 12]> year 4.43e+ 0 2.90e+ 0 1.53 0.266
## 3 grizzly <gropd_df [4 × 12]> year 6.1 e+ 0 5.16e+ 0 1.18 0.359
## 4 mid_teton <gropd_df [6 × 12]> year 1.21e+ 0 2.13e+ 0 0.569 0.600
## 5 n_teton <gropd_df [6 × 12]> year 1.90e+ 0 2.58e+ 0 0.736 0.503
## 6 paintbrush <gropd_df [5 × 12]> year 5.70e+ 0 3.70e+ 0 1.54 0.221
## 7 s_ak_basin <gropd_df [4 × 12]> year 5.71e- 1 2.05e+ 0 0.279 0.806
## 8 s_cascade <gropd_df [3 × 12]> year -1.09e-15 1.72e-16 -6.35 0.0994
## 9 s_teton <gropd_df [6 × 12]> year 5.54e- 1 4.53e+ 0 0.122 0.909
## 10 skillet <gropd_df [3 × 12]> year 5.79e+ 0 2.85e+ 0 2.03 0.291
## 11 windcave <gropd_df [3 × 12]> year 1.55e+ 0 8.25e+ 0 0.188 0.882
Non-significant increases at all sites except south cascade
Visualization:
summer_2.5_withps <- merge(summer_3plus, summer.2.5.lms)
summer.2.5 <- ggplot(summer_2.5_withps, aes(x = year, y = dd_2.5, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >2.5 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.2.5
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 row(s) containing missing values (geom_path).
Days over 5 C:
summer.5.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(dd_5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.5.lms
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 12]> year 2.50 1.44 1.73 0.333
## 2 delta <gropd_df [4 × 12]> year 0 0 NaN NaN
## 3 grizzly <gropd_df [4 × 12]> year 6.7 5.08 1.32 0.318
## 4 mid_teton <gropd_df [6 × 12]> year -0.0286 0.188 -0.152 0.887
## 5 n_teton <gropd_df [6 × 12]> year 2.09 2.33 0.896 0.421
## 6 paintbrush <gropd_df [5 × 12]> year 9.30 4.66 2.00 0.140
## 7 s_ak_basin <gropd_df [4 × 12]> year 0.357 1.27 0.282 0.805
## 8 s_cascade <gropd_df [3 × 12]> year -1.61 2.46 -0.656 0.630
## 9 s_teton <gropd_df [6 × 12]> year 1.05 4.27 0.245 0.818
## 10 skillet <gropd_df [3 × 12]> year 6.21 2.35 2.64 0.230
## 11 windcave <gropd_df [3 × 12]> year 1.68 0.729 2.31 0.260
Non-significant increases at all sites except south cascade, mid teton, and delta (0 all years)
Visualization:
summer_5_withps <- merge(summer_3plus, summer.5.lms)
summer.5 <- ggplot(summer_5_withps, aes(x = year, y = dd_5, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >5 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.5
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 row(s) containing missing values (geom_path).
Days over 10 C:
summer.10.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(dd_10~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.10.lms
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 12]> year 0 0 NaN NaN
## 2 delta <gropd_df [4 × 12]> year 0 0 NaN NaN
## 3 grizzly <gropd_df [4 × 12]> year 6.1 e+ 0 6.01 1.02e+ 0 0.417
## 4 mid_teton <gropd_df [6 × 12]> year 0 0 NaN NaN
## 5 n_teton <gropd_df [6 × 12]> year 2.01e+ 0 3.20 6.30e- 1 0.563
## 6 paintbrush <gropd_df [5 × 12]> year 5.30e+ 0 3.70 1.43e+ 0 0.247
## 7 s_ak_basin <gropd_df [4 × 12]> year -5.93e-17 0.327 -1.81e-16 1
## 8 s_cascade <gropd_df [3 × 12]> year 0 0 NaN NaN
## 9 s_teton <gropd_df [6 × 12]> year 1.57e+ 0 3.54 4.43e- 1 0.681
## 10 skillet <gropd_df [3 × 12]> year 0 0 NaN NaN
## 11 windcave <gropd_df [3 × 12]> year 0 0 NaN NaN
Nothing significant.
Visualization:
summer_10_withps <- merge(summer_3plus, summer.10.lms)
summer.10 <- ggplot(summer_10_withps, aes(x = year, y = dd_10, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >10 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.10
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 row(s) containing missing values (geom_path).
Make and save a table with stream names, years of data available for each, and stream source:
First, table set-up:
stream_info <- as.data.frame(summer_3plus) %>%
dplyr::select(full_name, year, source)%>% #dplyr::select full name, year, and source columns
group_by(full_name, source)%>% #group by full name and source
summarise(range = paste(min(year, na.rm = T), max(year, na.rm = T), sep = "-"))%>% #add date range column for each stream
arrange(source, full_name)%>% #sort by source then stream name
mutate(source = case_when(source == "glacier"~"Glacier", source == "snowmelt"~"Snow melt", source == "sub_ice"~"Subterranean ice"))
## `summarise()` has grouped output by 'full_name'. You can override using the
## `.groups` argument.
Convert to FlexTable and export:
tbl1 <- flextable(stream_info)%>%
theme_vanilla()%>%
align(part = "all", align = "left")%>%
set_header_labels(full_name = "Stream Name", source = "Hydrologic Source", range = "Date Range")%>%
set_table_properties(layout = "autofit")
tbl1
Stream Name | Hydrologic Source | Date Range |
|---|---|---|
Cloudveil | Glacier | 2019-2021 |
Delta | Glacier | 2017-2021 |
Middle Teton | Glacier | 2015-2021 |
Skillet | Glacier | 2018-2021 |
Grizzly | Snow melt | 2018-2021 |
North Teton | Snow melt | 2015-2021 |
Paintbrush | Snow melt | 2017-2021 |
South Teton | Snow melt | 2017-2021 |
S. AK Basin | Subterranean ice | 2016-2021 |
South Cascade | Subterranean ice | 2015-2021 |
Wind Cave | Subterranean ice | 2015-2020 |
save_as_docx(tbl1,path = "manuscript/tables/table1_streaminfo.docx", align = "center")
save_as_image(tbl1,path = "manuscript/tables/table1_streaminfo.png", res = 300)
## [1] "manuscript/tables/table1_streaminfo.png"
#Data set-up:
mean.long <- mean(summer_3plus$long)
mean.lat <- mean(summer_3plus$lat) #calculate mean lat/long of all sites
map_dat <- as.data.frame(summer_3plus)%>%
dplyr::select(full_name, source, lat, long)%>%
distinct()%>%
filter(!row_number(.)==10)
# Basemap: use get_map to pull map from google API:
basemap <- get_map(c(mean.long,mean.lat), #center map on mean lat and long of all sites
source = "stamen", zoom = 11, maptype = "terrain-background") #set map type and zoom
## ℹ <]8;;https://maps.googleapis.com/maps/api/staticmap?center=43.741346,-110.824899&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxxhttps://maps.googleapis.com/maps/api/staticmap?center=43.741346,-110.824899&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxx]8;;>
## ℹ Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.
bb <- attr(basemap, "bb") #extract max/min lat and long of the basemap for clipping
##GTNP boundary:
gtnp <- rgdal::readOGR("gis_data/gnp_gtnp")%>% #read in stored shapefile of GTNP and GNP polygons
spTransform(crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) #Update CRS to lat/long (comes in as UTM)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Users/GordonG/Library/R/x86_64/4.2/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Users/GordonG/Library/R/x86_64/4.2/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/GordonG/iCloud Drive (Archive)/Documents/Work/hotaling_data_analyist/tasr/teton_trends/gis_data/gnp_gtnp", layer: "gnp_gtnp"
## with 2 features
## It has 20 fields
## Warning: PROJ support is provided by the sf and terra packages among others
gtnp_clip <- crop(gtnp, extent(bb[1,2], bb[1,4], bb[1,1], bb[1,3]))%>% #Clip GTNP shapefile to the extent of the basemap
fortify() #fortify to convert from SpatialPolygon to dataframe for plotting
## Regions defined for each Polygons
Mapping:
site.map <- ggmap(basemap)+ #basemap: terrain background (no labels) from google maps platform console, centered on avg lat/long of all sites
geom_polygon(data = gtnp_clip, aes(x = long, y = lat, group = group), color = "black", linetype = "dashed", fill = "darkgreen", alpha = 0.25, size = 0.6)+ #Adding clipped polygon of GTNP
geom_point(data = map_dat, aes(x = long, y = lat, fill = source, group = NULL), color = "black", size = 3, pch = 21)+ #add sites, color code by source
labs(color = "Hydrologic Source", fill = "Hydrologic Source")+ #updating legend title
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snow melt", "Subterranean ice"))+
theme_void()+ #removing lat/long from x and y axes
theme(legend.title = element_text(size = 9, face = "bold"), #adjusting legend title/text/plot title size, face, and position
legend.text = element_text(size = 9),
legend.position = c(0.84, 0.12), #Moving legend into the plot area
legend.background = element_rect(fill = "white", color = "black", size = 0.5), #adding box behind legend
legend.margin = margin(t=5, r=5, b=5, l = 5, unit = "pt") #Increasing margins on legend box
)+
ggrepel::geom_label_repel(data=map_dat, aes(x=long, y = lat, label=full_name, fill = source), size =3.25, position = position_dodge(width = 1), fontface = "bold", show.legend = FALSE)+ #adding labels for each stream, color-coded by source
geom_text_npc(aes(npcx = 0.84, npcy = 0.65, label = "Grand Teton \n National Park"), size = 4.5, fontface = 4, color = "black", hjust = 0.5, alpha = 0.25)+ #add label for Grand Teton NP
geom_text_npc(aes(npcx = 0.2, npcy = 0.8, label = "Targhee \n National Forest"), size = 4.5, fontface = 4, color = "black", hjust = 0.5, alpha = 0.25) #add lable for Targhee NF
site.map
First, group by source and re-calculate max, mean, min, and degree days:
summer_source <- stream_summer %>%
group_by(source, year)%>% # group by source type and year
summarise(t_xbar = mean(summer_xbar),
t_max = mean(max_xbar),
t_min = mean(min_xbar), #calculate mean daily temp, mean daily min and max
dd2.5_xbar = mean(dd_2.5),
dd5_xbar = mean(dd_5),
dd10_xbar = mean(dd_10),#calculate mean no. days with max temp >2.5, 5, and 10 C.
n_sites = length(unique(site)))%>% #calculate no sites for each year
filter(!is.na(year), !is.na(source)
#, !n_sites < 2
)%>%
mutate(source = case_when(source == "glacier"~"Glacier", source == "snowmelt"~"Snow melt", source == "sub_ice"~"Subterranean ice"))
## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.
Visualization:
source.means <- ggplot(summer_source, aes(x = source, y = t_xbar, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "", y = "Mean daily summer (July-Aug) temp., C")
source.max <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "Source", y = "Mean daily max. summer (July-Aug) temp., C")
source.min <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "", y = "Mean daily min. summer (July-Aug) temp., C")
source.means|source.max|source.min
Snow melt streams look to be significantly warmer and have higher max and mins - no surprise there. Check significance of these differences:
max.aov <- aov(t_xbar~source, data = summer_source)
summary(max.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## source 2 55.74 27.868 22.95 2.73e-05 ***
## Residuals 15 18.21 1.214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(max.aov)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = t_xbar ~ source, data = summer_source)
##
## $source
## diff lwr upr p adj
## Snow melt-Glacier 3.68929231 2.036898 5.341687 0.0000981
## Subterranean ice-Glacier -0.08553797 -1.737933 1.566857 0.9900886
## Subterranean ice-Snow melt -3.77483028 -5.427225 -2.122436 0.0000769
Snowmelt significantly warmer than sub-ice or glacier, no difference between sub-ice and glacier fed.
First, visualization - re-organize data for easier plotting:
source_plot <- as.data.frame(summer_source) %>%
dplyr::select(source, year, t_xbar, t_max, t_min)%>%
pivot_longer(-c(source, year), names_to = "measure", values_to = "temp")
Plotting:
ylab <- expression(bold("Stream Temperature " (degree*C)))
source.summer <- ggplot(source_plot, aes(x = year, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Date", y = ylab, color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
source.summer
## `geom_smooth()` using formula 'y ~ x'
Test for significance:
Nest plot data by source and measurement type:
source_nested <- source_plot%>%
nest(data = -c(source, measure))
Modeling:
source.lms <- source_nested %>%
mutate(model = purrr::map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
source.lms
## # A tibble: 9 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> year 0.175 0.0225 7.77 0.00148
## 2 Glacier t_max <tibble> year 0.134 0.0514 2.61 0.0594
## 3 Glacier t_min <tibble> year 0.181 0.0283 6.38 0.00309
## 4 Snow melt t_xbar <tibble> year 0.135 0.412 0.328 0.760
## 5 Snow melt t_max <tibble> year 0.337 0.554 0.609 0.575
## 6 Snow melt t_min <tibble> year 0.0364 0.306 0.119 0.911
## 7 Subterranean ice t_xbar <tibble> year -0.0517 0.0986 -0.524 0.628
## 8 Subterranean ice t_max <tibble> year -0.000823 0.131 -0.00629 0.995
## 9 Subterranean ice t_min <tibble> year -0.0694 0.107 -0.647 0.553
Significant increases in min, mean, and max summer temp in glacier-fed streams, no change in snowmelt fed streams, non-significant increases in sub_ice streams.
Calculate degrees of warming per decade for each:
warming_amounts <- source.lms %>%
mutate(deg_decade = 10*estimate)%>%
dplyr::select(source, measure, deg_decade, p.value)%>%
rename(Source = source, Measure = measure, Degrees.Warming = deg_decade, P = p.value)%>%
mutate(Measure = case_when(Measure == "t_xbar"~"Mean", Measure == "t_min"~"Minimum", Measure == "t_max"~"Maximum"))
warming_amounts
## # A tibble: 9 × 4
## Source Measure Degrees.Warming P
## <chr> <chr> <dbl> <dbl>
## 1 Glacier Mean 1.75 0.00148
## 2 Glacier Maximum 1.34 0.0594
## 3 Glacier Minimum 1.81 0.00309
## 4 Snow melt Mean 1.35 0.760
## 5 Snow melt Maximum 3.37 0.575
## 6 Snow melt Minimum 0.364 0.911
## 7 Subterranean ice Mean -0.517 0.628
## 8 Subterranean ice Maximum -0.00823 0.995
## 9 Subterranean ice Minimum -0.694 0.553
Convert to FlexTable and export:
tbl2<-flextable(warming_amounts)%>% #convert to flextable
theme_vanilla()%>% #update theme to "vanilla"
set_table_properties(align = "center", layout = "autofit")%>% #center align the table on the screen, autofit columns
align(part = "all", align = "left")%>%#left align all text
set_header_labels(Source = "Hydrologic Source",
Measure = "Temperature Type",
Degrees.Warming = "Change per decade, C",
P = "P-value")
tbl2
Hydrologic Source | Temperature Type | Change per decade, C | P-value |
|---|---|---|---|
Glacier | Mean | 1.749728195 | 0.001477039 |
Glacier | Maximum | 1.340152287 | 0.059442199 |
Glacier | Minimum | 1.805660683 | 0.003094522 |
Snow melt | Mean | 1.350701233 | 0.759600210 |
Snow melt | Maximum | 3.371404238 | 0.575377043 |
Snow melt | Minimum | 0.364068864 | 0.911137685 |
Subterranean ice | Mean | -0.517086607 | 0.627706676 |
Subterranean ice | Maximum | -0.008229723 | 0.995279993 |
Subterranean ice | Minimum | -0.694147310 | 0.552638240 |
save_as_image(tbl2, path = "manuscript/tables/degrees_warming.png", res = 300) #save as .png @ 300dpi
## [1] "manuscript/tables/degrees_warming.png"
save_as_docx(tbl2, path ="manuscript/tables/degrees_warming.docx", align = "center") #save as .png @ 300dpi
Adding these P-values to the plot:
source_pvals <- merge(source.lms, summer_source)%>% #add model info to source dataframe
mutate(ypos = case_when(measure == "t_max" ~ 0.95, measure == "t_xbar" ~ 0.9, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
source.summer.pv<- source.summer+
geom_text_npc(data = source_pvals, aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.summer.pv
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_summer.pdf", source.summer.pv, device = "pdf", units = "in", height = 5, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
Re-organize data for easier plotting and modeling:
sourcedd_plot <- summer_source %>%
dplyr::select(source, year, dd2.5_xbar, dd5_xbar, dd10_xbar)%>%
pivot_longer(-c(source, year), names_to = "dd_type", values_to = "no_days")
Test for significance:
Nest plot data by source and measurement type:
sourcedd_nested <- sourcedd_plot%>%
nest(data = -c(source, dd_type))
Modeling:
sourcedd.lms <- sourcedd_nested %>%
mutate(model = purrr::map(data, ~lm(no_days~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
sourcedd.lms
## # A tibble: 9 × 8
## # Groups: source [3]
## source dd_type data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier dd2.5_xb… <tibble> year 2.43 1.30 1.86 0.136
## 2 Glacier dd5_xbar <tibble> year 1.79 1.30 1.38 0.241
## 3 Glacier dd10_xbar <tibble> year 0 0 NaN NaN
## 4 Snow melt dd2.5_xb… <tibble> year -2.95 4.01 -0.735 0.503
## 5 Snow melt dd5_xbar <tibble> year -2.08 3.80 -0.546 0.614
## 6 Snow melt dd10_xbar <tibble> year -1.36 3.03 -0.449 0.676
## 7 Subterranean ice dd2.5_xb… <tibble> year -0.518 3.32 -0.156 0.883
## 8 Subterranean ice dd5_xbar <tibble> year 1.00 1.04 0.963 0.390
## 9 Subterranean ice dd10_xbar <tibble> year 0.0714 0.169 0.423 0.694
No significant changes in degree days in any group.
Plotting with these P-values:
sourcedd_pvals <- merge(sourcedd.lms, sourcedd_plot)%>% #add model info to source dataframe
mutate(ypos = case_when(dd_type == "dd2.5_xbar" ~ 0.85, dd_type == "dd5_xbar" ~ 0.9, dd_type == "dd10_xbar" ~ 0.95)) #create column for adjusting P-value label Y position
sourcedd.summer<- ggplot(sourcedd_pvals, aes(x = year, y = no_days, color = dd_type))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","orange", "darkorange"), labels = c(">10 C", ">2.5 C", ">5 C"))+
labs(x = "Year", y = "Num. degree days", color = "Degree day threshold")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = dd_type, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
sourcedd.summer
## `geom_smooth()` using formula 'y ~ x'
First, combine glacier and snowmelt streams and recalculate mean, min, and max summer temperatures:
lumped_source <- stream_summer %>%
mutate(source_lumped = ifelse(source == "sub_ice", "sub_ice", "other"))%>% #add new col with lumped sources
filter(!is.na(year), !is.na(source))%>% #remove rows with no year or source
group_by(source_lumped, year)%>% #group by lumped source and year
summarise(t_xbar = mean(summer_xbar), #calculate mean, min, and max summer temps
max_xbar = mean(max_xbar),
min_xbar = mean(min_xbar))%>%
pivot_longer(!c(source_lumped, year), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
## `summarise()` has grouped output by 'source_lumped'. You can override using the
## `.groups` argument.
Nest data by lumped source and measurement type:
lumped_nested <- lumped_source%>%
nest(data = -c(source_lumped, measure))
lm for each source + year combo:
lumped.lms <- lumped_nested %>%
mutate(model = purrr::map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
lumped.lms
## # A tibble: 6 × 8
## # Groups: source_lumped [2]
## source_lumped measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 other t_xbar <tibble> year 0.220 0.290 0.760 0.490
## 2 other max_xbar <tibble> year 0.345 0.413 0.835 0.451
## 3 other min_xbar <tibble> year 0.146 0.204 0.716 0.513
## 4 sub_ice t_xbar <tibble> year -0.0517 0.0986 -0.524 0.628
## 5 sub_ice max_xbar <tibble> year -0.000823 0.131 -0.00629 0.995
## 6 sub_ice min_xbar <tibble> year -0.0694 0.107 -0.647 0.553
Add p-values to data and plot:
lumped_pvals <- merge(lumped_source, lumped.lms)%>% #add model info to source dataframe
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "max_xbar" ~ 0.95, measure == "min_xbar" ~ 0.85)) #create column for adjusting P-value label Y position
lumped.summer<- ggplot(lumped_pvals, aes(x = year, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source_lumped, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Year", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
lumped.summer
## `geom_smooth()` using formula 'y ~ x'
ggsave("Temperature/plots/lumped_summer.pdf", lumped.summer, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
First, calculate April 1 SWE and summer air temp from snotel data; add to mean temperature dataset:
#Calc means across both sites:
snotel_means <- snotel %>%
mutate(date1 = as.Date(date, format = "%m/%d/%y"), #recode date column as R date
airtemp_c = as.numeric(airtemp_c))%>% #recode airtemp as numeric
group_by(date1)%>% #group by date
summarise(swe_xbar = mean(swe_cm), #calculate mean SWE and airtemp across both stations
airtemp_xbar = mean(airtemp_c, na.rm = T))%>%
mutate(mo = month(date1), yr = year(date1)) #Add month and year cols
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `airtemp_c = as.numeric(airtemp_c)`.
## Caused by warning:
## ! NAs introduced by coercion
#Calculate max SWE and mean summer airtemp for each year:
snotel_summary <- snotel_means %>%
group_by(yr)%>% #group by year
mutate(swe_max = max(swe_xbar))%>% #extract max SWE and add to new col
ungroup()%>% #ungroup to get all columns back
filter(mo == 7 | mo == 8)%>% #extract summer months
group_by(yr)%>% #group by year
summarise(airtemp_summer = mean(airtemp_xbar), #calculate mean summer air temp
swe_max = unique(swe_max))%>% #keep max SWE
rename(year = yr)
#Merge with long-term temperature averages:
stream_snotel <- merge(summer_3plus, snotel_summary)
LMs - first, nest by site:
stream_snotel_nested <- stream_snotel %>%
nest(data = - site)
lm for each site:
swe.temp.lms <- stream_snotel_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "swe_max") #remove intercept term to just get info on year for each site
swe.temp.lms
## # A tibble: 11 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 windcave <tibble [3 × 14]> swe_max -0.00146 0.0127 -0.115 0.927
## 2 s_cascade <tibble [3 × 14]> swe_max -0.0165 0.00937 -1.77 0.328
## 3 mid_teton <tibble [6 × 14]> swe_max 0.000784 0.00325 0.241 0.821
## 4 n_teton <tibble [6 × 14]> swe_max -0.0561 0.0405 -1.39 0.238
## 5 s_ak_basin <tibble [4 × 14]> swe_max 0.00414 0.00603 0.687 0.563
## 6 s_teton <tibble [6 × 14]> swe_max -0.0522 0.0492 -1.06 0.348
## 7 delta <tibble [4 × 14]> swe_max -0.00178 0.00155 -1.15 0.370
## 8 paintbrush <tibble [5 × 14]> swe_max -0.0634 0.0543 -1.17 0.327
## 9 skillet <tibble [3 × 14]> swe_max -0.0305 0.00316 -9.67 0.0656
## 10 grizzly <tibble [4 × 14]> swe_max -0.0567 0.0840 -0.675 0.569
## 11 cloudveil <tibble [3 × 14]> swe_max -0.0166 0.0329 -0.504 0.703
Add P values and merge:
swe_summer_ps <- merge(swe.temp.lms, stream_snotel)
swe.summer.temp <- ggplot(swe_summer_ps, aes(x = swe_max, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Annual maxium SWE, cm", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'
lm for each site:
air.temp.lms <- stream_snotel_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "airtemp_summer") #remove intercept term to just get info on year for each site
air.temp.lms
## # A tibble: 11 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 windcave <tibble [3 × 14]> airtemp_su… 0.112 0.273 0.410 0.752
## 2 s_cascade <tibble [3 × 14]> airtemp_su… -0.251 0.0964 -2.60 0.233
## 3 mid_teton <tibble [6 × 14]> airtemp_su… 0.0643 0.0829 0.775 0.481
## 4 n_teton <tibble [6 × 14]> airtemp_su… -0.457 1.32 -0.346 0.747
## 5 s_ak_basin <tibble [4 × 14]> airtemp_su… 0.348 0.250 1.39 0.299
## 6 s_teton <tibble [6 × 14]> airtemp_su… 1.49 4.16 0.359 0.738
## 7 delta <tibble [4 × 14]> airtemp_su… -0.214 0.0826 -2.59 0.122
## 8 paintbrush <tibble [5 × 14]> airtemp_su… -0.0494 5.21 -0.00947 0.993
## 9 skillet <tibble [3 × 14]> airtemp_su… -0.346 2.58 -0.134 0.915
## 10 grizzly <tibble [4 × 14]> airtemp_su… 3.18 5.68 0.559 0.632
## 11 cloudveil <tibble [3 × 14]> airtemp_su… 1.81 1.04 1.74 0.332
Add P values and plot:
air_summer_ps <- merge(air.temp.lms, stream_snotel)
swe.summer.temp <- ggplot(air_summer_ps, aes(x = airtemp_summer, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'
First, add SNOTEL summary info to mean summer temps grouped by source:
source_snotel <- merge(summer_source, snotel_summary)
Re-organize data for analysis and plotting:
source_swe <- source_snotel %>%
dplyr::select(source, swe_max, t_xbar, t_min, t_max)%>% #dplyr::select temperature, source, and swe cols
pivot_longer(!c(source, swe_max), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
Nest data by source and measure:
sourceSwe_nested <- source_swe %>%
nest(data = -c(source, measure))
LMs for each source and measure:
source.swe.lms <- sourceSwe_nested %>%
mutate(model = purrr::map(data, ~lm(temp~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "swe_max") #remove intercept term to just get info on year for each site
source.swe.lms
## # A tibble: 9 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> swe_max -0.00523 0.00895 -0.585 0.590
## 2 Glacier t_min <tibble> swe_max -0.00349 0.00961 -0.363 0.735
## 3 Glacier t_max <tibble> swe_max -0.00689 0.00801 -0.861 0.438
## 4 Snow melt t_xbar <tibble> swe_max -0.0491 0.0355 -1.38 0.239
## 5 Snow melt t_min <tibble> swe_max -0.0417 0.0239 -1.74 0.156
## 6 Snow melt t_max <tibble> swe_max -0.0545 0.0532 -1.02 0.364
## 7 Subterranean ice t_xbar <tibble> swe_max 0.0192 0.00656 2.92 0.0431
## 8 Subterranean ice t_min <tibble> swe_max 0.0231 0.00566 4.07 0.0152
## 9 Subterranean ice t_max <tibble> swe_max 0.0106 0.0139 0.765 0.487
Add p-values and plot:
sourceSwe_pval<-merge(source_swe, source.swe.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
source.swe<- ggplot(sourceSwe_pval, aes(x = swe_max, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.swe
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_swe.pdf", source.swe, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
Re-organize data for analysis and plotting:
source_air <- source_snotel %>%
dplyr::select(source, airtemp_summer, t_xbar, t_min, t_max)%>% #dplyr::select temperature, source, and swe cols
pivot_longer(!c(source, airtemp_summer), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
Nest data by source and measure:
sourceAir_nested <- source_air %>%
nest(data = -c(source, measure))
LMs for each source and measure:
source.air.lms <- sourceAir_nested %>%
mutate(model = purrr::map(data, ~lm(temp~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "airtemp_summer") #remove intercept term to just get info on year for each site
source.air.lms
## # A tibble: 9 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> airtem… 0.249 0.221 1.13 0.322
## 2 Glacier t_min <tibble> airtem… 0.300 0.219 1.37 0.243
## 3 Glacier t_max <tibble> airtem… 0.124 0.229 0.542 0.616
## 4 Snow melt t_xbar <tibble> airtem… -0.690 1.12 -0.616 0.571
## 5 Snow melt t_min <tibble> airtem… -0.731 0.780 -0.937 0.402
## 6 Snow melt t_max <tibble> airtem… -0.437 1.61 -0.272 0.799
## 7 Subterranean ice t_xbar <tibble> airtem… 0.147 0.315 0.468 0.664
## 8 Subterranean ice t_min <tibble> airtem… 0.312 0.321 0.972 0.386
## 9 Subterranean ice t_max <tibble> airtem… -0.117 0.410 -0.284 0.790
Add p-values and plot:
sourceAir_pval<-merge(source_air, source.air.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
source.air<- ggplot(sourceAir_pval, aes(x = airtemp_summer, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Mean July-Aug. air temp., C", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.air
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_air.pdf", source.air, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
First, calculate richness, mean density, and total abundance for each site for each year; remove sites with only 1 year of data:
oneYr <- density_source%>%
group_by(site)%>% #group by site
summarise(nYrs = length(unique(Year)))%>% #count no. years
filter(nYrs <3) #extract sites w/ only 1-2 yrs data
invert_means <- density_source %>%
group_by(stream_code, full_name, site, Year)%>% #group by stream and year
summarise(richness = length(unique(taxa_lwr)), #calculate yearly richness, density, and total abundance
density_xbar = mean(density_xbar),
tot_abund = sum(abundance_total),
source = unique(source))%>%
filter(!site%in%oneYr$site) #drop sites with 1-2yrs of data
## `summarise()` has grouped output by 'stream_code', 'full_name', 'site'. You can
## override using the `.groups` argument.
First, nest data by site:
invert_nested <- invert_means %>%
nest(data = -site)
LM’s of Richness ~ year for each site:
richness.yr.lm <- invert_nested %>%
mutate(model = purrr::map(data, ~lm(richness~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
richness.yr.lm
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [6 × 7]> Year -1.34e+ 0 0.956 -1.41e+ 0 0.233
## 2 cloudveil <gropd_df [3 × 7]> Year 1.50e+ 0 0.866 1.73e+ 0 0.333
## 3 delta <gropd_df [5 × 7]> Year -1.10e+ 0 0.907 -1.21e+ 0 0.312
## 4 grizzly <gropd_df [5 × 7]> Year 1.00e- 1 0.790 1.27e- 1 0.907
## 5 mid_teton <gropd_df [7 × 7]> Year -8.93e- 1 0.394 -2.27e+ 0 0.0725
## 6 n_teton <gropd_df [7 × 7]> Year -1.54e+ 0 0.886 -1.73e+ 0 0.144
## 7 paintbrush <gropd_df [5 × 7]> Year -8.00e- 1 0.542 -1.48e+ 0 0.236
## 8 s_cascade <gropd_df [6 × 7]> Year -2.93e+ 0 0.852 -3.44e+ 0 0.0264
## 9 skillet <gropd_df [5 × 7]> Year -1.01e-13 1.25 -8.08e-14 1.00
## 10 s_teton <gropd_df [7 × 7]> Year -1.43e+ 0 0.818 -1.75e+ 0 0.141
## 11 windcave <gropd_df [7 × 7]> Year -7.86e- 1 0.340 -2.31e+ 0 0.0686
Add p-values and plot:
richness_withps <- merge(invert_means, richness.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
richness.yr <- ggplot(richness_withps, aes(x = Year, y = richness, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Taxonomic richness", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
richness.yr
## `geom_smooth()` using formula 'y ~ x'
Significant decline in richness at S Cascade, non-significant declines at most other sites
Save:
ggsave("invert_data/plots/richness_year.pdf", richness.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LM’s of Density ~ year for each site:
dens.yr.lm <- invert_nested %>%
mutate(model = purrr::map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
dens.yr.lm
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [6 × 7]> Year 12.9 12.2 1.06 0.351
## 2 cloudveil <gropd_df [3 × 7]> Year 64.6 16.9 3.82 0.163
## 3 delta <gropd_df [5 × 7]> Year 15.9 41.0 0.389 0.724
## 4 grizzly <gropd_df [5 × 7]> Year 6.06 44.9 0.135 0.901
## 5 mid_teton <gropd_df [7 × 7]> Year 42.9 39.5 1.09 0.327
## 6 n_teton <gropd_df [7 × 7]> Year -7.10 18.6 -0.382 0.718
## 7 paintbrush <gropd_df [5 × 7]> Year -12.9 203. -0.0637 0.953
## 8 s_cascade <gropd_df [6 × 7]> Year -2.25 16.2 -0.139 0.896
## 9 skillet <gropd_df [5 × 7]> Year 32.1 30.6 1.05 0.372
## 10 s_teton <gropd_df [7 × 7]> Year 3.68 2.05 1.79 0.133
## 11 windcave <gropd_df [7 × 7]> Year 28.6 19.6 1.46 0.205
Add p-values and plot:
density_withps <- merge(invert_means, dens.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
dens.yr <- ggplot(density_withps, aes(x = Year, y = density_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
dens.yr
## `geom_smooth()` using formula 'y ~ x'
No significant trends in average density at any site
Save:
ggsave("invert_data/plots/density_year.pdf", dens.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LM’s of Density ~ year for each site:
abund.yr.lm <- invert_nested %>%
mutate(model = purrr::map(data, ~lm(tot_abund~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
abund.yr.lm
## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [6 × 7]> Year 60.5 149. 0.406 0.705
## 2 cloudveil <gropd_df [3 × 7]> Year 333. 133. 2.51 0.242
## 3 delta <gropd_df [5 × 7]> Year -110. 383. -0.287 0.793
## 4 grizzly <gropd_df [5 × 7]> Year 107. 235. 0.455 0.680
## 5 mid_teton <gropd_df [7 × 7]> Year 4.89 190. 0.0258 0.980
## 6 n_teton <gropd_df [7 × 7]> Year -28.2 109. -0.259 0.806
## 7 paintbrush <gropd_df [5 × 7]> Year -383. 347. -1.10 0.350
## 8 s_cascade <gropd_df [6 × 7]> Year -193. 65.0 -2.97 0.0411
## 9 skillet <gropd_df [5 × 7]> Year 269. 281. 0.955 0.410
## 10 s_teton <gropd_df [7 × 7]> Year -137. 96.4 -1.42 0.215
## 11 windcave <gropd_df [7 × 7]> Year 50.6 119. 0.424 0.689
Add p-values and plot:
abund_withps <- merge(invert_means, abund.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
abund.yr <- ggplot(density_withps, aes(x = Year, y = tot_abund, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Total abundance", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
abund.yr
## `geom_smooth()` using formula 'y ~ x'
No significant trends in total abundance at any site.
Save:
ggsave("invert_data/plots/abundance_year.pdf", abund.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
First, group by source and calculate mean richness, density, and abundance:
source_inverts <- invert_means %>%
group_by(source, Year)%>%
summarise(richness_xbar = mean(richness),
density_xbar = mean(density_xbar),
abund_xbar = mean(tot_abund))
## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.
First, nest data by source:
source_inverts_nested <- source_inverts %>%
nest(data = -source)
LMs of richness ~ year for each source:
ricness.source.lm <- source_inverts_nested %>%
mutate(model = purrr::map(data, ~lm(richness_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
ricness.source.lm
## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [7 × 4]> Year -0.310 0.341 -0.909 0.405
## 2 snowmelt <tibble [7 × 4]> Year -1.69 0.895 -1.89 0.117
## 3 sub_ice <tibble [7 × 4]> Year -1.51 0.537 -2.80 0.0379
Add p-values and plot:
richnessSource_pv <- merge(source_inverts, ricness.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
richness.source <- ggplot(richnessSource_pv, aes(x = Year, y = richness_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean taxonomic richness", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
richness.source
## `geom_smooth()` using formula 'y ~ x'
Significant decline in richness at subterranean ice sites
Save:
ggsave("invert_data/plots/richness_source.pdf", richness.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LMs of density ~ year for each source:
dens.source.lm <- source_inverts_nested %>%
mutate(model = purrr::map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
dens.source.lm
## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [7 × 4]> Year -12.1 27.2 -0.444 0.676
## 2 snowmelt <tibble [7 × 4]> Year 0.528 30.4 0.0174 0.987
## 3 sub_ice <tibble [7 × 4]> Year 10.7 12.2 0.875 0.422
Add p-values and plot:
densSource_pv <- merge(source_inverts, dens.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
dens.source <- ggplot(densSource_pv, aes(x = Year, y = density_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
dens.source
## `geom_smooth()` using formula 'y ~ x'
No significant trends in average density in any source.
Save:
ggsave("invert_data/plots/density_source.pdf", dens.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LMs of abundance ~ year for each source:
abund.source.lm <- source_inverts_nested %>%
mutate(model = purrr::map(data, ~lm(abund_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
abund.source.lm
## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [7 × 4]> Year -103. 159. -0.649 0.545
## 2 snowmelt <tibble [7 × 4]> Year -110. 85.0 -1.29 0.254
## 3 sub_ice <tibble [7 × 4]> Year -47.4 101. -0.470 0.658
Add p-values and plot:
abundSource_pv <- merge(source_inverts, abund.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
abund.source <- ggplot(abundSource_pv, aes(x = Year, y = abund_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
abund.source
## `geom_smooth()` using formula 'y ~ x'
No significant trends in average abundance in any source.
Save:
ggsave("invert_data/plots/abundance_source.pdf", abund.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'